setwd("C:/Users/86139/Desktop/R作业/ADdata/")
library(openxlsx)
a <- read.table("ADdata1.txt")
b <- read.csv('ADdata2.csv')
c <- read.xlsx("ADdata3.xlsx")
d <- read.table("ADdata4.txt")

colnames(b)[1] <- "ProteinID"
colnames(c)[1] <- "ProteinID"
merged_data1 <- merge(a, d, by.x=0, by.y=0, all=TRUE) 
colnames(merged_data1)[1] <- "ProteinID"
merged_data2 <- merge(b, c, by="ProteinID", all=TRUE)
merged_data <- merge(merged_data1,merged_data2,by="ProteinID")

ADdata <- merged_data[,-1]
#e <- log2(ADdata)
ADdata1 <- as.matrix(ADdata)
f <- log2(ADdata1)
g <- pmax(f,0)
#修改列名 删除列名后两位数字
colnames(g) <- substring(colnames(g), 1, nchar(colnames(g))-2)
Group <- colnames(g)
groups <- c("asym", "ad", "ctl")
ProteinID <- merged_data[,1]
# 创建一个存储结果的数据框
result <- data.frame(ProteinID = character(0), P_value = numeric(0))
# 循环处理每个蛋白
for (i in 1:nrow(g)) {
  # 提取第 i 行蛋白数据
  pro_data <- g[i,]
  
  # 判断各组数据是否有至少3个样本是有数值的
  pro_ctl <- Group == "ctl"
  ctl <- pro_data[pro_ctl]
  pro_ad <- Group == "ad"
  ad <- pro_data[pro_ad]
  pro_asym <- Group == "asym"
  asym <- pro_data[pro_asym]
  
  # 计算各组非0的个数
  ctl_num <- sum(!is.na(ctl))
  ad_num <- sum(!is.na(ad))
  asym_num <- sum(!is.na(asym))
  
  # 对于判定结果进行方差分析或者直接输出NA
  if (asym_num >= 3 & ad_num >= 3 & ctl_num >= 3) 
  {
    anova_result <- oneway.test(pro_data ~ Group)
    p_value <- anova_result$p.value
  } 
  else 
  {
    p_value <- NA
  }
  result <- rbind(result, data.frame(Protein = ProteinID[i], P_value = p_value))
  result$P_value[is.nan(result$P_value)] <- NA
}


#GO富集#####
library(ggplot2)
library(org.Hs.eg.db)
library(dplyr)
library(clusterProfiler)
library(forcats)
# 导入数据
load("C:/Users/86139/Desktop/R作业/volcano.RData") 
data <- prostat
geneID <- data$ID
FC <- data$FC
P <- data$P
pro_genes <- geneID[P < 0.05]
enrich_result <- enrichGO(gene = pro_genes, 
                          OrgDb = org.Hs.eg.db, 
                          keyType = "UNIPROT", 
                          ont = "all", 
                          pvalueCutoff = 0.05, 
                          qvalueCutoff = 0.2, 
                          readable = TRUE)
# 创建气泡图
result_BP<-enrich_result%>%filter(ONTOLOGY=='BP')
result_CC<-enrich_result%>%filter(ONTOLOGY=='CC')
result_MF<-enrich_result%>%filter(ONTOLOGY=='MF')
BP<-result_BP[1:10,]
CC<-result_CC[1:10,]
MF<-result_MF[1:10,]
all<-rbind(BP,CC,MF)

ggplot(all, aes(x = -log10(pvalue), y = fct_reorder(Description, -log10(pvalue), .fun = median, .na_rm = TRUE), group = ONTOLOGY)) +
  geom_point(aes(size = Count, fill = p.adjust), shape = 21, color = 'black', na.rm = TRUE) +
  facet_grid(ONTOLOGY ~ ., scale = 'free_y', space = 'free_y') +
  scale_fill_gradient(low = 'blue', high = 'red') +
  labs(title = 'GO Enrichment', y = 'GO term', x = '-log10(pvalue)') +
  guides(fill = guide_colorbar(reverse = TRUE, order = 1)) +
  theme_bw()